{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multitask GP Regression\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This notebook demonstrates how to perform standard (Kronecker) multitask regression.\n",
    "\n",
    "This differs from the [hadamard multitask example](./Hadamard_Multitask_GP_Regression.ipynb) in one key way:\n",
    "- Here, we assume that we want to learn **all tasks per input**. (The kernel that we learn is expressed as a Kronecker product of an input kernel and a task kernel).\n",
    "- In the other notebook, we assume that we want to learn one tasks per input.  For each input, we specify the task of the input that we care about. (The kernel in that notebook is the Hadamard product of an input kernel and a task kernel).\n",
    "\n",
    "Multitask regression, first introduced in [this paper](https://papers.nips.cc/paper/3189-multi-task-gaussian-process-prediction.pdf) learns similarities in the outputs simultaneously. It's useful when you are performing regression on multiple functions that share the same inputs, especially if they have similarities (such as being sinusodial). \n",
    "\n",
    "Given inputs $x$ and $x'$, and tasks $i$ and $j$, the covariance between two datapoints and two tasks is given by\n",
    "\n",
    "\\begin{equation*}\n",
    "  k([x, i], [x', j]) = k_\\text{inputs}(x, x') * k_\\text{tasks}(i, j)\n",
    "\\end{equation*}\n",
    "\n",
    "where $k_\\text{inputs}$ is a standard kernel (e.g. RBF) that operates on the inputs.\n",
    "$k_\\text{task}$ is a special kernel - the `IndexKernel` - which is a lookup table containing inter-task covariance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up training data\n",
    "\n",
    "In the next cell, we set up the training data for this example. We'll be using 100 regularly spaced points on [0,1] which we evaluate the function on and add Gaussian noise to get the training labels.\n",
    "\n",
    "We'll have two functions - a sine function (y1) and a cosine function (y2).\n",
    "\n",
    "For MTGPs, our `train_targets` will actually have two dimensions: with the second dimension corresponding to the different tasks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_x = torch.linspace(0, 1, 100)\n",
    "\n",
    "train_y = torch.stack([\n",
    "    torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * 0.2,\n",
    "    torch.cos(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * 0.2,\n",
    "], -1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the model\n",
    "\n",
    "The model should be somewhat similar to the `ExactGP` model in the [simple regression example](../01_Simple_GP_Regression/Simple_GP_Regression.ipynb).\n",
    "\n",
    "The differences:\n",
    "\n",
    "1. We're going to wrap ConstantMean with a `MultitaskMean`. This makes sure we have a mean function for each task.\n",
    "2. Rather than just using a RBFKernel, we're using that in conjunction with a `MultitaskKernel`. This gives us the covariance function described in the introduction.\n",
    "3. We're using a `MultitaskMultivariateNormal` and `MultitaskGaussianLikelihood`. This allows us to deal with the predictions/outputs in a nice way. For example, when we call MultitaskMultivariateNormal.mean, we get a `n x num_tasks` matrix back.\n",
    "\n",
    "You may also notice that we don't use a ScaleKernel, since the IndexKernel will do some scaling for us. (This way we're not overparameterizing the kernel.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MultitaskGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(MultitaskGPModel, self).__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = gpytorch.means.MultitaskMean(\n",
    "            gpytorch.means.ConstantMean(), num_tasks=2\n",
    "        )\n",
    "        self.covar_module = gpytorch.kernels.MultitaskKernel(\n",
    "            gpytorch.kernels.RBFKernel(), num_tasks=2, rank=1\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultitaskMultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "    \n",
    "likelihood = gpytorch.likelihoods.MultitaskGaussianLikelihood(num_tasks=2)\n",
    "model = MultitaskGPModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the model hyperparameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/50 - Loss: 47.568\n",
      "Iter 2/50 - Loss: 42.590\n",
      "Iter 3/50 - Loss: 37.327\n",
      "Iter 4/50 - Loss: 32.383\n",
      "Iter 5/50 - Loss: 27.693\n",
      "Iter 6/50 - Loss: 22.967\n",
      "Iter 7/50 - Loss: 18.709\n",
      "Iter 8/50 - Loss: 13.625\n",
      "Iter 9/50 - Loss: 9.454\n",
      "Iter 10/50 - Loss: 3.937\n",
      "Iter 11/50 - Loss: -0.266\n",
      "Iter 12/50 - Loss: -5.492\n",
      "Iter 13/50 - Loss: -9.174\n",
      "Iter 14/50 - Loss: -14.201\n",
      "Iter 15/50 - Loss: -17.646\n",
      "Iter 16/50 - Loss: -23.065\n",
      "Iter 17/50 - Loss: -27.227\n",
      "Iter 18/50 - Loss: -31.771\n",
      "Iter 19/50 - Loss: -35.461\n",
      "Iter 20/50 - Loss: -40.396\n",
      "Iter 21/50 - Loss: -43.209\n",
      "Iter 22/50 - Loss: -48.011\n",
      "Iter 23/50 - Loss: -52.596\n",
      "Iter 24/50 - Loss: -55.427\n",
      "Iter 25/50 - Loss: -58.277\n",
      "Iter 26/50 - Loss: -62.170\n",
      "Iter 27/50 - Loss: -66.251\n",
      "Iter 28/50 - Loss: -68.859\n",
      "Iter 29/50 - Loss: -71.799\n",
      "Iter 30/50 - Loss: -74.687\n",
      "Iter 31/50 - Loss: -77.924\n",
      "Iter 32/50 - Loss: -80.209\n",
      "Iter 33/50 - Loss: -82.885\n",
      "Iter 34/50 - Loss: -85.627\n",
      "Iter 35/50 - Loss: -87.761\n",
      "Iter 36/50 - Loss: -88.781\n",
      "Iter 37/50 - Loss: -88.784\n",
      "Iter 38/50 - Loss: -90.362\n",
      "Iter 39/50 - Loss: -92.546\n",
      "Iter 40/50 - Loss: -92.249\n",
      "Iter 41/50 - Loss: -93.311\n",
      "Iter 42/50 - Loss: -92.987\n",
      "Iter 43/50 - Loss: -93.307\n",
      "Iter 44/50 - Loss: -93.322\n",
      "Iter 45/50 - Loss: -92.269\n",
      "Iter 46/50 - Loss: -91.461\n",
      "Iter 47/50 - Loss: -90.908\n",
      "Iter 48/50 - Loss: -92.142\n",
      "Iter 49/50 - Loss: -93.466\n",
      "Iter 50/50 - Loss: -90.492\n"
     ]
    }
   ],
   "source": [
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},  # Includes GaussianLikelihood parameters\n",
    "], lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "n_iter = 50\n",
    "for i in range(n_iter):\n",
    "    optimizer.zero_grad()\n",
    "    output = model(train_x)\n",
    "    loss = -mll(output, train_y)\n",
    "    loss.backward()\n",
    "    print('Iter %d/%d - Loss: %.3f' % (i + 1, n_iter, loss.item()))\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make predictions with the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAd8AAADOCAYAAAB//QjSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzsnXlYVVXXwH/nXrj3Ml8QEByDzFLTRKnMJgcsNdM0TPO1bNL63uaM0tTKJsqy8i0bFMvKLBLnBlMhs5wSxFlzAGcmgct85/P9cQe5CIiCcNH9ex4e4Jx99tlnWGftvfbaa0myLCMQCAQCgaDxUDR1AwQCgUAguNwQylcgEAgEgkZGKF+BQCAQCBoZoXwFAoFAIGhkhPIVCAQCgaCREcpXIBAIBIJG5rJXvpIkaSVJmiBJUg9JkmIkSZpQaV+MJEmLGrk970mS9FIN+2IkSSqs0saXJEn6UpIkbTVl11yE9sXaf/eQJCnN3t7ISvsjHfesLvevcjurObZB2l/5ntqfd2xD1CtwP4Q8n3f7hDw3EZe98gUWAT/JsrxNluW1QIHjwdr/b2wSa9phb8+cKpu3ybL8uCzLumrK6mhAJEmKAbbZ698GZACJsixnVDpvhizLIyu1oVYqt7OaYxuq/c576rhPlT8wgksKIc91RMhz03JZK19JknrAmQdo/zsJmFypWKS91xZrf1kdvcSYSj9ae481xt7rjpEkaY39mJfsv9Ps5WIkSfqmcnl7nS/Z6485R7O/BB6v9L/Wfrzz3NVcZ+Xe6EuSJL1n/7tqu12uq5pzD6gsmDXd0+p6uPbrj616znMday/3kuNZ2bc52jrhHNuqvaf2Z1z5HgouAYQ8C3luTlzWyheIxtbbOwvpjNmnQJbltfYH/J592yhw9uYysAn3Wvv/Pe2/I2VZTpJleYb92AL7R6EAyK5c3v6COXrqtfYuHcIi2Uw6WiBDsvX6ImVZngO8XM0xlXudSZV2ubS7muuqiraabVXPtY0qPVy7wIys9CGsfM7ajo20l0tytE2yjWJSHT11u5BWt+1c97TZ9ZQF50TIs5DnZsPlrnwzqP6haSv1niu/PAV2wYgHBkiSdBjbC9wDCLL35r60l91Wpc5F9pc2qJryA6jho1EDjt5yjN28liHL8hypyjxRHajajqrXdV7U0LuOtLfVUV9196omqrsnAzjzTDLs/9e0rbZ7WnCOcwuaH0KehTw3Gy5r5WvvRQVVfskl2+T9jErFKr+0OntPNUaW5Zex9fRigDX2+hzzJtXxE7aXtqCa8luxvcRgE+ZztXsO4HQysJtnqnXqqIbKH6eq7ah6XedLj2q2pQLjOTPKqMu9qo1tnLmGSGz3rrpt53VPBc0fIc9CnpsTHk3dADegPzBBkiRHzzbI/sI6cJiBtJzp2V0vSZJjf5Isyxn2+QgA7L97SJLUw/5SIsuyTpKkAvv/2yqXl2V5RqV5kB7YeqtzqjpdVKGy+SUSW08x0t7eWGyC0EOSpEj7B2arvSerBWIkSdJWOq+jTpfrquaczvZUauso+/0ZAOgc2x2mM2ymwCD7vvdkWX65yr0qqFReW8PfMfa/tVWO7yHL8gx7PTVtq+meNqjzisBtEPIs5LlZIMkiq5GgjtiFPeNcThruzqVyHQJBfbhU5KC5XsdlbXYWnB92s151pqhmg8Mk2dwEVSBoaIQ8Ny1i5CsQCAQCQSMjRr4CgUAgEDQy9Xa4io6OdnjRDUhNTT1rTZpAIGgeCFkWCBqPeo18o6Oje2AT1LVAj+jo6Ga30FkgEAhZFggamwaZ842OjtYC76WmptYY4mvSpEliclkgqCPvvvuudO5SDY+QZYGgYalJlhtqnW80dVhnNX369Fr35+bmEhoa2kBNahhEm+qGaFPdqEubXnvttUZqTbU0iCxD873/jYm7tQdEm+pKfWW5QRyu7KYqbXR0dLNM7SQQCGwIWRYIGof6zvm+Fx0d7cg6oeMSCfslEFxuCFkWCBqX+pqdvwQi7V6S2tTU1Kq5KQWXAGazmezsbAwGAzX5CFitVoqLixu5ZbXj7m2SJAm1Wk1YWBgeHk0e6VXI8mWAkOWGo76yXC+JT01NzeBMMO2mSFQtaASys7Px8fGhdevWVIoV64LJZMLT07ORW1Y77t4mWZbR6XRkZ2fTpk2bJm2XkOXLAyHLDUd9ZVkE2RCcE4PBgFarrVFYBReGJElotVoMBkNTN0VwmSBk+eJwIbIslK/gnMiyLIT1IiFJUo3mP4GgoRGyfPE4X1kWylfQoGRlZRETE0N2dnZTN0UgENQDIcsXF6F8BQ1KfHw8Gzdu5J133rngOtLT01myZAnJyckkJCSQkWGbilyyZAlTpkxpqKbWiE6nY/DgwdW2q1OnTiQnJ5OcnMzMmTPR6ZpdGlGBoE4IWb64COUraBCCg4PRaDTMmTMHq9XKnDlz0Gg0aLXa86pHp9ORkJDAiBEj6N+/P4899hhTp04FoF+/fhej6Weh1WqJiIg4a3tUVBQRERH079+f/v37M3HiRMaMGVNtHTqdjpkzZ17spgoEDY6QZVculiw3+foGwaXB7t27mTp1KitWrKCiogIvLy+GDRvGu+++e171JCUlERUV5bItMDCQ9PR0IiIiSE9PJz09nZSUFB599FHS0tIICgoiJSWF4cOHk5KSQlBQEN27dyctLY2lS5cSERHB1VdfzapVq1i4cCFPPvkkEydOBHApHxQURFJSEhEREWRmZtapvVqtFp1OR0FBASkpKRQVFTF8+HAyMzNJS0sjPT2dgIAA5767776bjh07ntc9EQgaEyHLjSPLQvkKGoSwsDD8/f0xGAxoNBoMBgP+/v6EhYU16HmioqKcAj1v3jwACgoKiI2NZerUqTz66KMEBQUxb948nnvuOdatW8fbb78NwKpVqwAYMWIEkZGRjBkzxqW8Tqdj4sSJREZGsmTJkjq3qaCggMjISIKCgpwfiUcffZSUlBRnWx37li1bxksvvdSQt0QgaFCELDeOLAuzs6DByM3NZfz48axfv57x48eTk5Nz3nXExsaSkpLisi0zM/OsHrSD4cOH8+ijj/L+++9jNBrp2bMnUVFRTiGtbCrr168fM2fOpGfPns5tVcsHBdkCOxUWFtapvTqdjsjISJKTk50fkMpkZGTUuE8gcFeELF98WRYjX0GDkZiY6Px71qxZF1SHVqslLi6OhIQEp8no008/dSnjMFVNnDiRmTNn0q9fP2JjY+nevbvT1BQUFEReXh6ZmZnodDq0Wi2xsbFMmTLFKcRvvfWWS/mJEyeSlJREz549yczMJD093eVDkZGRQWZmJsnJyQBs377d2bbMzEwiIiIoKipylnNsr7yvunoFAndDyPLFl+UGSSlYFyZNmiSLrEYNQ2O36dChQ3To0KHWMu4egcZdqK5NVe/va6+91mQpBetCXWQZhOzUBSHLdaO5tOl8ZFmYnQUCgUAgaGSE8hUIBAKBoJERylcgEAgEgkZGKF+BQCAQCBoZoXwFAoFAIGhkhPIVCAQCgaCREcpX4FZkZGRw0003kZ6eXus2gUDg3ghZrh0RZENQZzQadS17a9t3Nnp99UmnIyMjnQvzZ8+eDUBRUREREREiMIVA0IDULM9ClhsDoXwFbkdAQECN+zIyMlwCqANnBUGfOXMmb7/9NomJiecdDF4gEDQcQpZrpl5m5+joaG10dHSs/ee9hmqUwD3R6w01/pSUlNa6v+rPuRgxYgRLliw5K3zb1KlTiYiIICIignnz5hEZGekMR7d06VL69++PVqslKiqKoqKii3k7LimELF9+CFluWuo753sfEJSampoEEB0dPaH+TRIIoH///iQlJVW7r3IA9eoCnQcGBjZGEy81hCwLLgpClqunXmbn1NTUOZX+jQS+rK18bm5urfXpdLr6NOeiINoEVqsVk8lUaxmLxdIg59q+fTtz586lW7dudO/eHR8fH7Zu3UpGRgZbt27l9ddfJzExkYiICAIDAzl06BDt2rWjoKCAQ4cO8d1335GRkeEMxL5161anScsdqO4+Wa3Wc8rGxaahZRmE7NQFIcuXryw3yJxvdHR0JFCQmpqaUVu5ugQQd6eg5w4u9zYVFxfXKah5QwQ+v/766/nxxx8BXHJlbt682fl35QTW119/vfPvkSNHAvDAAw8AsHLlSrcLxg5n3yeFQuE271hDyvL5lGtM3K1NQpbPlHdwOchyQy01ik1NTX28geoSCARNh5BlgaARqLfyjY6Ojk1NTZ1h/zum/k0SCARNgZBlgaDxqJfZ2S6g70VHR0+2b3q5/k0SCASNjZBlgaBxqa/D1VrgygZqi0AgaCKELAsEjYsILykQCAQCQSMjlK9AIBAIBI2MCC8pOC8++eNwtdstFitKZd36ck/3FdZNgcAdqE6ehSw3DmLkK3A70tPTSUhIIDk5mSVLljBlypQ6H5uQkEB6evpZx+h0OgYPHtzQTRUIBLUgZLlmxMhX4FbodDref/99Fi5c6NyWlpZW52MzMzN57LHHiIiIcNmn1WrP2iYQCC4eQpZrRyhfgVuRlJREv379XLbFxcWRkZHB9u3bCQgIoKioiICAAGfGk6SkJN5++23S0tLIzMwkOTmZmTNnsnz5cnQ6HUlJSURERJCZmQmcnU2lcvYUR10AM2fOJCIigoCAACIiIlyOiYyMbPR7IxA0J4Qs144wOwvcHq1Wy9SpU+nXr58zSHvljCeO+Lg9e/YkIiKC/v37O3vGU6ZMcR7n2FY1m0p1dSUkJNCvXz/69etHUVHRWccIBILzR8jyGYTyFbgVsbGxpKSkuGxLTk6utmxdM54EBQUBUFhY6NxWOZtKdXWlpKQQERGBVqtlxIgR1R4jEAhqRshy7Qizs+C8qMm70WQyNUjgc61WS1xcHAkJCURERFBUVES/fv2IiIhwmpxiY2NJT08nMzPTOTeUnp5OQUGBy7bMzEwmTpxIUlISPXv2dJZ76623nHU5hLlqXXFxcS5lqh5TOS+pQNBcqU6ehSyfmwqjBQVyve6PJMv1q6CuTJo0SZ4+fXqtZXJzc90u44hoExw6dIgOHTrUWqahBLYhaS5tqnp/X3vtNd59912psdtWV+oiyyBkpy4IWa4b7tQmk8VKYbmJALWEWqVy2Xc+sizMzgKBoMHJKdajN1mbuhkCQYNiscoUVZhoiEGrUL5VyMrKIiYmhuzs7HOW2blz5znLXgpIktQgL5vgbGRZRpLcdpB7wRSWm1h9oJBjJ05eFjLSXBCyfOHIskyx3oTFWv39O19ZFsq3CvHx8WzcuJF33nnnnGXGjRt3zrKXAmq1Gp1OJ4S2gZFlGZ1Oh1qtbuqmXBROl5n4b9yrl4WMNBeELNcdo9HIv//+S3l5Of/++y+FpRUYzdVbcy5EloXDlR2tVoter3f+P2fOHObMmYNaraaoqKjaMvv27XMpq9FonO7tWVlZPPDAAyxYsICwsLCzzneu/e5EWFgY2dnZFBQU1Ci0VqsVhcK9+nLu3iZJklCr1W7//C+EG65ug9FgcP5fnYwIGh8hy+fGYrGQn5+PUqmkrKyMnJxczGYTpwt0ePsFOMsVekp4KJXAhcmyUL529u3bx6RJk1ixYgUVFRV4eXkxbNgwXnzxxbPKLF++3EUJazQa7rnnHt59913ntsoj6P/9739nne9c+90JDw8P2rRpU2sZd3NkAdGmpuSXv1J5/oUX+Xfrn5gMelRqDbfGDOLNd+KxWGWUikvP1N4cELJ8bnx8fLBYLDUXkCSmLfiDMTe2J6JN+AWfx726N01IeHg4/v7+GAwGNBoNBoMBhdqbbLMP6w+eZtWeHLbkyJw2KjEYDEj2XphCYfu/yOJBZrmKgIAANBoNc+bMwWq1Onv7Wq0WsI2ea9svEFwKhISGofb2wWw04KFSYzIaKEfFhlMW5v59hOT9eRSWG5u6mQKBE8e3uSbF66FS06JVOyRg7cLP632+y175mi1WSvVmCsqMHDuZxX1jH2LGN8voPWQUOw8eI/VECTtOFHEwt5TjhRXk5OTSa/AoIrteT8v2HYjoGk2vwaM4fjKb1KOFvPTV70T1uQuVWgOAxsuLUaNGs3//fsA2eh41ahReXl62/RoNISEhrF+/vsnugUDQ0OTlZrN3UzI9+g/l6Y8W0mvwKApyTvFZ3Djy83LYm1XM91tO8NueHPJKDOeuUCC4SGRlZdGvfwxvJCwhqs9deNq/3VUxGw3knzqGLEts+kVHpw7f1WvQdEmbnc81r1qqN7N8RxYF9h54zNPvAWAAhv3fVADKysooLshjQfyLjJ38AeOmzar1nP5BIai9fTDZe/wGvZ4T5bBHp8TT13jWCFuv16PX60lISHB787NAUFfm/G8m+tJiVGoNrSKvYcRT01jy6Rts/vUnfv3qQwpyTjHsicnMjotn7OQPuCayHd1a+3NlsA+KSibp5uQbIWieTH39TTZt3IBB8wWHdmzBbDQgSQpk2Yp/i1C8fP0pKy6ktDAYGAc8CLQG4OefH7ng817Syre2edWCMiMrdmRRYjBXe6xD4Y549g02LvuGI3vSWLvwc4Y/+Somg4RRr8BskvD2s6DSuDoulOoK6DV4FL0Gj2Tzr4soyj/N3qxi9mYV42ctZVHSYiRJqtbBSzikCJozNqdEAG9AZtMviWz6JdGlTFryCgA+fupeJEli7cLP8X/qVU7pKvBVe9C9bQDXtQ5AoZCalW+EoHlR1YHW8V4CPD87ic2/LmLTL0kU5/cCngUqJ4k4SNdu2+nYcegFn7/eyjc6OjoGeDk1NXVAfetqKGryXJYkiczMTCTvQBb8sYOEN19g7OQP8A8KOauONd9/TuZuMzPHrwBuBf6PTb90YdMvQVS11nuqrfj4W9CGmAhpU0D20d7cMmwIYVf4MOKpaS5lv/nsQ4p0hfS+8x60GgUpv/9KRUUFGo0GPz8/fvnll4twRwSCc9MQsrxv3z4efDCRv/56GchGkvYSFKZD7Z3JqcPzgH0u5WX5jIL28FQRvyKdvw/lc1uX9piMwltacHEwmq18umw9s9+dTvq6s7+5Hz35AArFEwSGFlOY62XfWgok0vaabYS0PkmQUk9oaBMq39TU1LXR0dEv17eehqQ6z+XWrVtz+PBhnp/8Gn0ensTP3852jmZHPPUqAGYTTBk+CatlEPABUL1XoNLTgsZLRukhU16ixGRQoMtToMvz5Mheb2A6X7wMvoFmuvQq4drexcyf3gWLucxZx4ZVS51/q9RqYX4WNDkNIcvh4eEoFKHYPlRhyHIY+VmOvROB08DfQAqwAjiKh0pN15tjGPJYnLOeSV//zi8J77NnUwoG/ZnVB5VXFAgEF4Ku3MQvu7MplH1Qe/sAOM3M4INC+SwKRRxmk5bCXFB75WCoeA+l5wKs5tO0ufI+Rjz1DkOu8q5XOxrV7Jybm1vr/obq0SqVSjw9PZ2j34qKCg4dOgTA0h++YekP3zjL2nrdKUiK/6JQPInVsqJSTSeBZGAHCuW/WC3pQA5Rfe9hyATbN0qWwahXMOPhMVjM7YGeQDRwA6WFbdnyWyBbfgtE5VVIQMhvFJ2eisW0Cw+VGm+/AEoKT1e7HlKtVnPw4MGLep8aEtGmuuGObboQapNlL6+fGBqbjOqKYWxPOUhhVgC6vHAkxe2YjcHAPfaf/wHbMRuXYzTmo1R7U1Zm66Aq1d5InmqMBj0enir0ej0KpQcKheKc35HacLf7727tgebdppycHJ566ilmz55d7dKkEzoDf2YUYbAHy9CdziX6znspLSxj/z/dgDislmCsFmjdoZRbRmSx/Y8n8AsKJHrAR6SuWUphXg5lZWUUFZnIVSkv+JoaVfnWZZ1WQ63lKikpYcKECUTdOoCpL71AfvYJZKsVD5UaP20LSnT5mI3hSIrpyNZRyFY1Fit4+Z4gKGwzJw+9i9JzNxaTASQJq+XMvG7q6sWkrl7sNJPhC698M4ef577P7k0fYzLo8VBpuLLbY4S2ncjB9DCyj2goqBgKDEWSkjEbZ9H5Rm8G/Of/7MclYzLo0Wi8uOceWw+/tnvhbuvwQLSprrhjm86X2q5h2bJl7M8uYVlqJjf27+J0tLrxjvsoyPbEYu5N9tGrKC/uhSx3B7qzZwOUFZZzw52FdLu1GLWXjKG02MV3YteRbPbqFPRsp6W0MO+CHbHc7f67W3ug+bbprbfeYuvWrcyZM8fFgmi1ymzKLCD9pAkPtRce9kBUD786m7SUAJZ86gnY6vfVHqBF+Nc8OfMBJAmi+37irOfKLrYMSMUFeTz1xFMkLfrpgh0BLymHq8qekYmJiWw9UsjmzALad+rO6VPH8PBUYTEZ8VC1xGx8AXgK2aoGrNhMYLOoKE3h5CGbGeKZjxfx14qFlOTn4u3r71SQnmoN1/bu72Imc3g5O9Y1WkwGglrmMXRCOZDBwe3lfPXqbqyW0Vit/YH+7PjzX67trUHldeY4g0FPscUT38DgJrmHAkFD8db9t2A2nVnLu/lXm+OVh+dPxK9Ix2Q8zuEdPqQle7LjL3+O7PXhyF5vln3ekl6Dixj2xGy0ITaHSIfvxK6TRew5VcwfX7k6Ygmv6Mubmvx8NBoNx7Ly+H1vLtnFepdjDm33ZsWcMLIybUuLWkVWMPiRXDr2MCNJD9R6vrULPyctdWu9HAEvKeXr8Ix8++23uefJV9mbVQxAxu40ADrd2I+Swvs4sncYEGQ/6lvgDeCws56ovncx5LE4/INCGDL+ZXx8fFj8yXTMRgOeKjVmo4E2oS24M/oaygxmdBUmiipM6HKz8A0MZvSL8ezesJaSwtNOr+mglq2wmFcSPeAvWrZ7mz8Xt6Ck8GoSpoK3/yt0uekqBvynB5t/XcSxk1l8t/k4XVv707OdFq96mDYEgqbi2c+WkfL97Bo7rZ4qmWuuL2XvljeQrb/gH/QkxQVDMOpvYf2SFvy9PIiovkX0jT1Ny/Y2JT55aJSLQnd8ZB2hB4VX9OXJn3/+yZAhQygpKUGv1+Pl5cVdQ4Yy9plX+DH1JAbzmcAZRac9WDm3JTvW20JFakONDByXR1SfIs4VwbKm9+9CHAEbwts51vYrOjY1NTWpvvVdCFV7PXPnzmXu3LlVSl3Jrr9fwea5DFd0LkTt9ToZu+dgMtiOlRQKkGU03r74B4WgUipoEaihS/tgUqQKJkyYwKOPPsq8efPIzs4mur1tgbWt1/0IPa9uz7/bNlC4529ej/+AYwXlTOh3DbLVSuZuWyu2rv4K+AqlRwCDH/mXP34Kprz4OnZvvA5v/0IGP9wFH38LZquVP3cc5PH7X+SjL74iJuoql/WPAkFD09Cy7BcY7GINMhsNTtmCsz9kxQXvAe8BUcBLWC0jSVurZVtyAD1jirhjbC6T5692maZxYLXa5vCEV/Tlybx588jLywNsySP0ej0nyyUOlnoANsVrMcNfy1qwdmEwhgolnmor/Ued5rZ78/FU1Z5oIjxAQ5tAL95amEzi7Hed7199HAEbwts5CWgSpeugqnezo4d9+72PsG7R1+zc0AGr+S3ABw9VAcOfPMX1AxQs+TQb87YzC6q73hxDi6AQTGUFDO/eivAADfmn8wgNDWDp4kXO882a5Rpoo0OHDlgsFv7++28Akr6fT9L382ts75mRdT69BhXyx6Jg1i9pwT+rAtm90Y+7HskleoCOtQs/J2N3GvHvvE3RK+8Q0ymEQG9VjfUKBPXhYshy1TXvJYWnnfscinTXxrWYKy0r8lDtpevNC7h5aASpa67in98DSV2jZfs6f24eGoTSI9RF8VblzjvvpLCwkOzsbGGCvsSpOvACbOF/JYXLu3bysJpFH7Xi5GHbsqFrexfT9769/JzwNNff+QGe1Sw3VSkVXBPuR9dW/gT52L670e20/LM0mO1GA2q1GoPBgL+//wW9Z5eE2Tk8PByVlw96vd6lhx0Y2pnMvW9gNUfZS35P99tTuOGOicDZHwaFvohffkhwjjCzsrIYPXo0iYmJ1d7c6h68g2HDhiHLMmvWrKGiogIAhUKBLMt4+/o5e/9evlYGP5xLdIyOJbPDObzDh0Uft2LRx3uBDVQOVKD0VNG1ew/mfvaJ0/lAzHUJ3Jlx02ZxdUs/booMQje8LzklRnJLDBSUGVFIoai9fbCYjEgKBbLViqRQYDEZ0Xj70v4aP9pfk83t9+bz+7ehbP8zgD8XB4M0gzYdu+Dtt5xD6RuxWl1j8W7ZsoXi4mJhgr4MqGng5Zg2NBklkn8I5o9FLbBaFAQE64l9Jpdrri9lyacfn7XcFECpkOjWOoDo9lo0nq5Tfh5KBZYyHY+NH8/w4cNZvmzZBeeqviSU7/GCcnYcPOZUpOuXfEv6ulPs+6c9RXleeHiWMPChPeSf+tmlN+QIFalUSIyb/QnXhPm51BsfH88///xToxBXffBwRsG2bNkSWZZtCRoUCqxWK8OHDyc4OJhTWVncflUwW44UojfZPhyhbY08Hn+U7ev8WfpZCBWlMcAu4Gk8VEl0vbk/CoWSbSkriXtjBku/n4+XSikiAAncFg+FggGdQp1y5afxoG3QmbWRVqvM2k8rGP3Aw+zdv5+cnFx8tC0IbRNBQc5JPosbx9jJHxDcKoT/TDrJ7ffm8/XreooLOnPiwET8AgdjtY4E9ric12FuFiboS5fKgw5/f/+zBl7+QSGcOKThx/dbkXNMg82p9hOu7pnKN2/+6DLd4QzyolKz7J+D3HhFEH6amlVjYqLNcTA3N5d+ffte8DU0a+WblZXFiFFjuPu5d3lg6hlTsO70HejLHkZf5kXbjhU8MCWLwFB/YNpZdQR4eXJn51Ba+p8Jpl2b51xlIa4cp7mqgnX0hsaPH+8yT1zZZH11S19Sj+nYcaIIi1VGkiCqbzH/bnuTtLWDgWHAd5iNw0j/YwJQCMC6lYsI9D9jBq+tjQJBUxDo7cmwa4O4qkqHtjIKhcTSpDPvscUqc7SgnK17M3h69GBKCk87RyWu88P/AWZSUtgJSKdHv10c3TeS/OxM28J7O45Uny+88AIxMTEsWLDA7fLUCi6Mys61+zJPuFgwiwsKSElswe/fhWK1SMC/wCPARv75/UwdDmWtUmvoN/AuPv5gBld86u4rAAAgAElEQVS0bbwlVs1S+Tp6PV5B4WzfuhnFvJkU5Jzi2P5dWMxvAC/ZS37J8QPPMuMx2bYeF/BTe9CtTQAhviqCfdXVehLXlNu3ukn13NzcWhWsg+q2qT2V3HxlC65t5c8/RwoZfsNVmE2Oua85wEPALCAWiEbpMRaLeQMeKjXXRN+CyWjkwLYNyFariAAkcCta+muQ9Of3eVEqJHp0aO3S8XVOuXh4EtXnLrujy/d4qJLRBidw+tRdbEvpgY//GpCHIin2Oc3XBoMBD403CQkJTuvQ1KlTG/pSBY1Idc61AB6eKkY8NY0+sW/ww/ut+W2+zcLS5aY9ZO65E6P+NGZ7303t7YuhvBSz0YBarcFkMtKuZQuuaNu6Ua+lWSpfh4OTA1tAbBUwH7gfMAFP4KleyLW9z4StC/JRMaxbOL61mBTAdUR7rkl1hwkCqlewdSHAy5MBnULZun0Xz0yMY/Mfv9uXZvzIVVEKjux9mfLijljMKSiUkzAbP8YvMJiMXanIVisKpbJeE/8CgTtQkw+FQqFgybo03ot/u9I6+hyuivqO2Ge7kPhhKwpzIpEU6QS2/AIPz0/wCwwiY9dWFsz/ylmPsA41fxwDo+UrVqCvMse7Y70/SbPC0Zcr8VQXYDLcz+mThykvPulSh6G89MzfBj1KpZKcnJzGvpTmlc+35mTHAcDv2BRvMTAID9X3Lvb/8AANsVGtzql4HThGtMuXL2f8+PGN8nA6XdmOa9qGYjYZ8VSpMRn07N38FeXFXbGF4lNhtXwI/MimX34m55htbbLVYsFqtZKQkEBWVhYxMTEX7AQgEDQVjlzXSqWrNWrMmDF0a+3Hvk1ruGXQCGd+4JLCfK7sVs4Ln2Vw48BCZKsnBVlPow3exthJ3zD1uxSX/KxeXl7cc889ztzaguZHeHg4JUYr+ooKPDxVttgL6kBWL+jGgvg26MuVwBJMhg7Aauc3sjo0Gg2jR4/m8OHDLoOoxqJZKd/0nbvpNeDuKsmOw7AFau8DnMJXew89Y3xdBDQi2Idh14Wj9lTWWTklJiYya9YsOnfuzKxZsxrt4eTm5jJh/HhWJ6+j98ARBAS3xFOtAJ5F6TEGhbIMuA+Nzy48VFcD4KnW0LPfELbs2OfigCUQNCccFieLxYJSqUSSJDp37kxJSQnx8fEU6XRc0zqI2AE3E/v0q06HSY23ldhns3hk+jF8/M0c2ObLh/+NJPd4e5d1xnq9AV9fX2EdaqaYLVaS9+exaeNGADr36sN1tz3PjvVvsOW3QDw8rQx66CDdb/8KD1W5y7GODp1jzl+pVGI0GpvUWthszM6lejMbs2Tw9KqU7DgchXI9VksksBcYRNebb3K6jY985lVuigyie5sAJMm2fMjdvYMrK/lVSQsY8eB4kpcutJvafuS621py/OA75J+KADag8BiL2fg729b9yvVdf3YeWzlBQ1FRURNciUBw/uTm5vL44487fSjmzp3L3r17nfsdAXTUajWzV++msPyM12qnG0p5fnYG37/XmszdPsx5pT0hrQdx4yAlN90Vy+ZfF3HouLAINUfKjBauDAx0STO5868A4HXAj5A2BsZOPkGrSBOLP3FdvqZUKrFYLHTu3JnQ0FByc3MJDQ3l6quvblILYbNQvmUGM18lb+fz155DpfGi1+BRdLpxHN/HR2GoaIm3/xEenJLBjvW9nUuJgn3V3NEplBa+tsXRdfVgbkqqrtlVeSjwNpdy130PcNVtw+xBCnYzPn4fs58zUlLYC9n6C+07f4Pa62uX+NOOF27kyJFNfVkCQZ2p6kMxadKkGp0fW4SEsjmzkF0nbasFAAKCzTz+7lHWfB9Cyo/B5B5/iODW9xLY8iQjnppGWVkZP/21iznTnxdr45sJ2UV6lu/OZ9LXv9uDsqzHbHwDeA6ALjflcH9cAWov2ztQqiug95DRWAtPUlRw2kXRNoV5uSbcXvkaTBZW7Mxm6VefcGRPGr0G38dtI97gy8ntMVSoaNuxgsfe0uPtF8mV3WxLibq1CeCWK1ugrBSO8Xw8mJuK6kblP/1ke1l2niii9ZWdkGWZsrIypn7nz+/fnSYlMZgjex7m5qFDMZuedEb+ccyLL1iwgAULFrhVJ0MgqCuVnR81Gs1ZjoW3dmhBVJsA0o7p2JNVjMUqo1TCwAfzuKJzOQvfa8PezX588lwED716HJ+gMj7/+AM2b9jA5FemcPzYUaGE3YTqAgYdzitj9d5c9GYr/kEhSIrWmI2/ATcDRq7o8j3jpt2IVCny7vjXP2HodeGEB2iqPY+74NbK12SxEtwiyMXUsOmXzWz6RQOoaN+5nEffOIaXjy2uq4TEzR2CiGqrPauucwlxU1KXUXm3NgFovT1ZtSeXMkChhEEP5RLazsCij8LZsKIFfoHP0f12qCg7xcFttsg/VTsZIiKWoLlR3XK+yvhqPLi9YzA922lZf+g0h/NsOYGviS7jmVkZfPNmW7KPaHh/QijwPGCbnvlh4fcAXHnllc48woKmo+rgI/24jg2HCpCxjWiP7PNi51/vAC3w1VZwRed4YBOSdKOzjvAADbd2CKalv7ppLuI8cFuHK4tVZtWeXCZ9/Xslj8VWQArQlrYdixj/1lGn4lVIEgM6h1areB04hHj9+vWN5sFcFxxenl5etrijXl5ejB49+iyvzHZB3sT2aIW/+kyfqWe/IibEH8Pb30xJYS9yjs3H268Lsmy1O5no8fXzQ5ZlYmJimDZtmnDIEjQrHM6P3bp1q9X50VfjweBrw+h9ZQunj0dwKxNPf5TJdbcVAf7AcuAFl+MsFgsajQattuZvh+Di4VjFMmfOHKxWq3PgcWvndk7Fm7YmmM/j2mM2taDdNTpe+OwY46aNZviTU/ksbhxe5hLuuS6c2B6tm4XiBTdTvpn5epZs2MN1vW5j5oqtHMkvc+bJNRn8gWQgEp+Aw0x4J8tp4/dUKhh6XThXt/Sttf66CnFjcz6j8iAfFdcHm0mY/DDFBbYsHpHXlvPwa2l4qo6Slalh14Y36HrLi06P7x0HjtGhQwf+/vtvFixY4PKCiw+O4FKjZzstw64Lx8sel1elkfnPpJPc+WAutk/eTCAB8ARsqwVuHTiMXXv21lSl4CJSdfCh0miI6nsXk+evxmKGxZ+EsfLLCKwWBfAJ4RFP4hdom1ZL+fELjuzZxs6VX7mELm0OuJXyTTtZyicfzuDAjq389t1s5/ai00a8/bcA1+Dtd4y2HePQ2Ee8SoXEXde2pG2gVxO1umE4n1H5l7P/x8GdqfyzJMG5bVvKR5iMPfAL3I/ZGMKBtHjKi3uwdfUStv75ezVro6l2dC0QXAq0DfSiTxsP5tk7qZIEMfefplWH15EURuBRbLEBgjAbDZiVGv4+ZaXUYHapR6ybv/hUHnx4qtQY9XoO7fiH3OMlvD6qgs2/BgF6bBH/nmHLbwuJG9SFuEFd2LDyR2S5eQ4m3GbOt+q855mwcgG07XiK8mJvWoSX4xMwnpHPvQLY5njv6BTa7Ho81VGXSFlV79Hvixfw++IFLmVKCrsDC9CXx5IwrT13T9jO0b1xLvlPHckf3GXOWyC4GPzvwxkc2JnK3l/nc+u4lzFZrEx49y4KTpzgs7gAzKa+BLXMpH3nqZQU7uN0qYGkbacY2i3MmULO3ZcmXiqcysqhz9D7uW7ACBa+9zI5xzTMnXIdVktbPFUFtO30Osf3J2Iy2CwV19/ch2BflTNrnDs60J4Ltxn57tu3jx63D3QG0PBUa+h++91Edj3Mkb3eaENMtO8Ux/F/17B24ecA3N4xmA6htZuaLyUc5hmN5kzEnmHDhjFw8N2onPdNonufb7lh4AksZolln3UiP3uUc200wJC7hzFhwgS3mfO+XCg3nm19EDQ8VecQF38/n+fu6MQrw3oA0PZqPS/PyyfsCj0FOf4cTP+Qfvd9CcDJU6e4tU8//P0Dqp2HbE4jq+bC/uwS7nx2But/WcRHT95LzrFOwEaslrbAVkzGrmTs+ASTwZG5yMi1HdrRsmVLt3SgrStuo3zDw8PRePs6o9GYDAZOHn6eg+ktgHx0ed3YlvIZsmzLbxs3qAu3dmnX1M1uVKqbG27ZsiVtW4dhNhlRq21ZOrx8vBn5bDFDHrOZyo7tf4CwK1by7P8WcdNdozlWUMpb7810mznvy4ESvZnUo4VN3YzLgpocGP/dv5+u4T4AaEPM/PeDI1wVVUqpzoPPX76CPZt9WbvwczJ2p9Ht1oEMGDLcWYdGoyEkJIT169c32XVdClQ241cYLfy2J4c1+3IxmC1M+no1rSK/ApYBfsD3wG14qPLxbxHKjXcOZ8WqtUyYYJuWc1cH2rriNmZngBJdvjM1VOKHSk4d7oun2srYyUfZvu5Kdm86gsmgR63RMPyee5qViaGhyM3NZezYsTz99NMuyy4cSzG++HIu2w8cBeD2ewvwDzKT+GFrsjLvYt3iIka9cA0enrAk/RQjolrVmrdS0HD8cSAP72oyaAkanpocGCUJZrzwCPGffsWOAgVePlYemX6MxZ+0InWNlvnTW2OLEy+TmrzcWZ9arUGv16PX60lISBDmZzsXsmzRYcZ/5bU3uOXBlykzmikuyOO7tyfj32IJpzJutxWUpoD8DpJCgcUkc12v21m68Gt81R706339WfVeaFKbpsStvrwPT5mJRalmw8pATh0OR6GQeeCVE3S6QcO+f2wxWj1VakxNHJOzKUlMTHSGR6vuhfts9ifoTRZW7swmu1hPVN9ifLQWvn2zDdvXBVBWpOTBqScAE4vSTjKkWxihfs3DNb+5si+rhKP55XQKrzm3raBhqW5tcHx8PP/88w8/zp3Fq++8z8+7cjhZcIrTp8Zx46Bv2fJbJ+BLoA1Kz7dBtmK1WDAY3DsyXlNxPvPhVf1VFn7zFQu/+QoPTxXdbx/Hkb0zgDAkRQVXRX2GxZxMqa4DvtoWRLRvS9pfayktPI3vJfTNr7fZOTo6OjY6OjomOjr6pXOXPjd7Nvuy/AvbDY597hStO2TyWdw4CnOz6DNsDOv+/LNZmhgaE42nknu6h9PO7ojWMaqMJ2YcxVdr5mC6L7MntuaT558nKzuLJemnnEEJavLsFB6fF06pwcxfh/Kbuhl1oqFluSmpvKzwm2++Yfny5cyZMwdZlpkzZw5tggN4dmBX/kz8kqN701AophPZ9UvADEzDYvoMi1mmR78hLpmRhPm55nW5tc2H79u3j3tH3uf0TfFQ2Tr8ZtNVpK6dDNwCHEe29iZj5ys88e5XTJqzksXLfyYyLAidTnfJxSaol/KNjo6OBUhNTV0L6KKjo2PqU9/Rf735/t02yFaJO8bmcv2AItYu/Jwje9IIDW/DT/O/pGdUd7dao+uueCoVDOkaRkf72uc2HfQ8OTOTFuFGso/4cWz/Z6ycuwSTxcpvu3PYdkx3VkYkh9IVgTmqpy6dkj8PnMZgdn9Hq4aWZXeiujlghUKBwWBg/YofnH4kGbueAO4ByoFHgOWkJSeTvu4XW35te9CavLw8EhISajnjpU1dgwJVRu3fgly9ApPdp8dsNACDkaTNQASwBQ/VbUT1bc3k+avRenkyaWh3ul8R4tJpupSc3uprdr4ecGjBDKAHsLamwrm5uTVWdPSokoTpkZgMCrr3zSP5x7asXnAmY8lfK38gRPsDarWagwcP1rPZdccdTUvn06buLcBUDruzy5j5xC2YTQHAr0A029e9wvZ1g4FtLsc4TGsO/v77b5ft1T2D5n6fLoRp06axYcMGpk6delbHxGSxsuVYCQfyKpzbinRWdAr5orapHjSYLDtwl3dCqVTi6emJXq9HpVKh1+sZPnw4FouF33//Hb3e5kWr0nhz73MD2bB8Mhk7pgKDQVqHr3YcZUX7XcLc1iYL54O73KPKnKtNle+nWm3rkHh4eKBQKKp9L06XmVj9byF5OdkgSXbF+xzwAbKsBH4AHsFiMqD07ElkWBC3tlfx999/89ZbbzmfkUajYeDAgUydOrVO79/Fpr7Prr7Kt2oXpEVthUNDQ6vdrtPBI494UlasoGOPUka/kEefe39g7pTxGCrKbB68ldZx1VTPxaKxz1cX6tqmrKwsPp70OJNmfMYr89ewcu4Mdm0chNn4PXAHCuXfDH18G0f3xrusBa6J0aNH1/gMmvN9Oh+qzl9VTV6RU2wgZW8Ox07msyD+RYY9MZnlX8TzwrS3GP/yayQmJrqjv0KDyPKFlrvYlJSUMGHCBIYPH87SpUvJzs4mNDQUo9GIRmNzqDIbDRxK+5vQ1jIZO24BVoEcjcmQzH/fP8rGlW86ZaQhv0fuco8qc642Oe6nY079yJEjjBkz5iznqxOFFfx1PBul2otHp8+mMPc0X75cQX72nbYC0usEhc3n3qc/Zc/GtXgairj/5o4AtAkPcz4jtVqN0WgkNDSUa6+99qJd9/lSn2dX3zlfHRBUzzrw84NBg6y0iijngVdOoPSALb/9RFlRAWZj813H5Q44TMkr5n/KyFu74OXji8WUj9LzXuA7rBYvVnzZm+KCQc5lXgChrds7zUpgC8whSZJ4BtRsdtu7dx9px3QkbTuJrsLknDL5YcZLHNmTxozJT/PPP/+4q/m+QWTZXXHMAXfu3Nk5bZWbm4skSWcF99n8609I0iEmxO/Fxz8DQ3kY89/oicnYzSkjer0eH1+/y1YWqobqbd++/VlTU4fzylixMwujxRaNsKxIyY/v97Qr3goUyrFIvMHVPXpx3Q238F3C5/y+cqnLeRyOc8uXL7/kfH3qO/LdypkecySw5kIqUSphxgwL4bceZNr9N2A2GV326/W2/LSX0o2/2NSUKUmhUHDb0Pvpcce9bPplEYd2eHD65P0c3vkkba9uwb1Pq9jy2yL2p/6FXq9HoVBgtVoZPnw4wcHBl63TVeVlFdUtZfHQePPnSQunS/OZPDTK5R3OOXYY0HI847/AfHf1mG0QWW5OJCYmkpWVdVaq0R633UG/B5/HPyiESV8b+e6dUg6k+bJn86t07BnI4IevZPOvi0g/cJRSgxlftetn9HLKHFbTd0alVvPOsnRnYoTso2q+fr0tBdkqPFSn6XT9R8SMGcHmXy2UF+UzIqqVM6pYZRy+Pbm5ufTp06dRrqmxqNfINzU1NQmIdDhn2J01Lhi1l5XJ81cT1ecup1ecY1Rx+PBh4WR1HtQ0OsvIyGDpd3OJ7hHFgP88gV/ghwwcdxBJkjn+72g2/dyXYU9Mo/WVneg1eBRxny1h7EOPYjabL2tHt6rOaI4e+ZqUPxgUO5btB45xutQ2J+h4h21WBCXwf8BBbPNcM9FoNG4XV7uhZbm5UF1HqlO7UK7vFAmAxtvKI68f4/oBhchWDQe2vcjRfTcx4qlp3D/5I5LSTpJf6jpYqPquXMpU953pM/geXv7qd2RkigvyeP/xeXzyXHsKslW0uaqCSfMKeHDqf2gVeQ1Dxz2Jl7UcY0lBE19J41Pvdb6pqakzGqIhDpxZjExGYW6uB+fKlHTPdeGM/PA1juxJI6z9NB6Y8j4LZ7Rmy6pACvM8eWDyJ87kFeGRL3LXtS1rPV9OTk61cz7NndpyLR88nsPPu7LpMewR9sa/yMnD+1n+RTxjJ3+A2tsHs/EW4COgq/3oP1Ao4jC66Tr1hpbl5kJ1a4IHdArBbJU5nFeK0gNGPp+FNtTMmu9DWPJpOIW5ngwcl0uJwcxXa7ez4uNJ7ExPw2A42ynLzSwcDUrl74zaPndukGxTV7NfHIex4jlyj80EoPvtRYx87hQqjW007Kv2YPOi+fyzedNlGTvbbYJsZGVl8emkx/jP5A+xlBUxoZbk2YK6UVMS8hqTWChvwSfgDw6k+fJZ3BU8Mv0Y2hAzJouVFTuzue2qYLq29q/2XLNmzbokA9Dv27fvLLPksGHDeCLuVRann8JksbrM7eYez+CXeYkc2vFfwBatR6E8jod6Ki3DNtM2LJiuXceLd9qNqCmpyZ2dQ1m5y8rxgnIkCe4Ym4c2xMTi/4Xzx0/B5By1UlYyguDwFmzbsom7ho/Ex1M661251CPx5ebmMnbcw7TrPZTkZT9QUniaN8feDfI8YCRgBV5h+5/vsXujivgV6Uwe1sPu9WzjcuioVMVtlG98fDyZe9JZu/BzFv74ozNFYHMMG+Yu1PRRqapQ1BovOt/UjyGPxWHUZ/LVq+3IytTwyXMRPPTqcdperccqy6w7kEd+mZHbOrRAobAlK69tZHgpCFF1FgS9pCY1FyYNva7K3O4x4EW2pbwK+AJlDBxXxm0jSvFUvUyncD+6Bbmnd6vgbJQKicFdWrJ8RxbZxbZ3/IY7dQQEm1jwThv2bgkFPuDo3mGAzM9LfnIeeylZ7c41h/3G/+aRsj8Ps9XK1tVLMJvaAxuBa4EiYAzwK1F972LIY3G08FGzfdce3nx1ymXXUalMkydWqBwtxbHY/arwwEtmIbU7UlWhmIwGIluFEBAUSnArE099lElk1zKKCzz57KUr2L7uzGh318kilu/MosKeoae6TEvuNp9ZH7KysliyZAljxozh59XJxAwfw6GjJ5GRq8zt3gZsB2YAvmhD/uaZWbvoP/o0niq3XdsrOAcqDwVDu4UR5q9xbps/vSP68uuAI0AvYAvQBYCA4Jb0GxJL8h/r6uSdm5OT4/bR42qaw7ZYZf46lM/qvTmYrbYpqkEPrwNSsSnefcCNSIpVSJKExtuXqyPacm9UOB3at6l1WuxyoMmVb9UJe43m0vp4uyuVM4KMGTOGdatW0llrU6jeflbGv32UGwcWYjYq+P69Nqz6NgS7fHGisIIfU09wSqc/59xycyc+Pp7CwkKKTQp2lAcy4LHJjJtmsyL4B4WgUIZjNn4J/Al0Bg4Ad9Lphk9p27F6E72geaH2VBLboxW3dwxG7aFk8vzVBLcqAW4ENgNXAJuAYXS+sQ+DnpzOXkMQ4596gfz8/FoVa+XpGnejtjCSunITi7edYvtxm3XLaoFV34awcs5NgD+wGFvH5F+63hxDr8GjsJbpGHZdOGpPW4KR5p6VqL40udm58sfbU6XCaLy0Pt7uSmWTtLe3N4WFhSz7+hNGPvMa24/r8PCEe5/JomV7AyvntiT5hxCyMjWMnngSL18rpQYzS7afoldEUI2Zltyd2sxpVc3pKxK/ZUXit3h42uasik7n8fnLGyjMmQn4I0kGvP2/ILTdYsLataak8HQjX43gYiJJEt1aB3Br5/YuiRagL/AVcD+wjMM7v8FqBV2FiWcmv8mmDa5+EI53buvWrW7vnFWTv8O9457g5tv7MmbSB/gHhVBc4MFbY08hy50BC/AKNgsQSJIC2WK1rQluE+BSf03TYpcLTa584UwPqO/Au/lj1cpm8/Fu7tQ0X+upUvPO8m1IEtx6TwEhbQwsfK8Nezf78fFT7fDye4RHpv8f/kEhbMrI5z9TPqZbEFzRJtwpRM1hrWNtWVl279nLhGde4K+1q2xxfdUaru3dnyGPxXEqQ83cKeGU6l4B4KqoUkY8lUVwqwHAgCa4EkFjsX+/TSEtWrQIq9WKh6cV/+DJKBQ68k89Qe7xcUwasgxZfhAoAVwV64MPPsjGjRsZM2YMJpOJ5cuXo9fr3XLO82xPZgP5Rg8+/NSW83jtws+5tvd7LHwvHFnuiIdnATJjsJh+d8rLuGenct+tXTAU5xMTc69bfw8amyY3O8OZaCm9e3a7rNeSNjY1rQVe/+efzJv8MMUFeQBcE13Go2+motIcoCDHm5OH5vHD+zsBKC7I4+VHRvLN+n/Zm1XsrLsuax2bKlvSubKy6MpNbMySKZdVzohGtnSWwbz1wDY+erI9pbqOQBYwioPpfsx8wn1C3gkuHg6FBLY8vxaziat79OLlhNt55I1jaHwsyPI9qL324eF5AwCeag2SQoFer3e+cwsWLCAxMdEZs9hdp2uOnsiiz9D7eXLmQpBgzZIFbPolEVn2YNMv1zF3SlvKilXAWrShg7GaVzvlpU1oC8bf0Z0gH9Vltfa5rriF8hU0DTXN13799Vcc2JnK1qVnMrfMntgbo7478C3gw6HtTxM3aBNv/mcIR/aksebHuSTvz8PPP6DO6caaSiBr6nSk7tjDH//m8f0/x8kp0VOqK6DX4FE89eFCOnT/iG1/fIBsfR5QoFB+BlyDp3oFUX3vYvL81S7n8FN7cPtVwdx+VbDIl3yJ4bDU/fXXeh5+5DGMJYWArZP6zKxMwiP0GCpaYzb9iUL5JCaDHm9/Ldf27u9MTejl5UXr1q2JjY112znPXSeLufPZGQx6/BXCI69m6ncpNgdDzyhsTmaTARl4FbiT0ye3IMsystXK/Q8+gkJfRIugwPNOP3i54BZmZ0HTUXkt8E033eSSzWhV0gJWJS2oVLoCGIdtGcFHwKPALcjy/aSuXkzq6sUoPTyJ6nsXezYlY6zBnNbUy5Oqdjr0BgO5pSaG3juKofYkCGMnf8C4abPIOaZi2edhHNpuS83YtmMFAcFvsGfTe3ioVJiNRjTevvgHhQAQ4OVJz3ZargnzQ2lfjtWtTQAFZUaK9WawlF706xNcXCpb5j6b/Qlmi5W/DuWz+1QxIa2NPP1RJu8/fpjCnCFYLZ+i9h5Ime4B8k4ccYkN3W/AQN5+4zVCQ0Pdas7TKsv8eeA0O08WuWz31YZQkBOL2TQa0AAZBLSIo7z0V0wGK55qDTf2uZO5n35IRNvWANWG73Q383pTIUa+lzmVA6QfOnTorBHhsGHDGHzX3c5wnwCSYi5wPUqP/cDVwGYUylfofvvdvPLNGtRePpgMZz4yeHq5mNMuJB9oQ1DZzJ2VncOI+8cx9Yskeg26j+1bN5NZKQnCqm++ZfkXLfnwv1dyaLsv3n5mYp85xbhp/5Cx6wt69B/K0x8tpNfgUZQU5uOhUND7yhaMvaEtXVr5OxWvgyAfFVe08L6o1ydoGjyUCkVjedoAABdiSURBVPpeHcLALi1RKRV4qmVemR+JQvkgUIqhfAiwm5xjHewjQwu9Bo9i+4GjbDtRitFsbepLcFJQZmTV/kJ2niyiuCCPz+LGUVyQR85RFZ/HXcHRfQ8BGrrcdIQb7nwThceWSlMzRjq1C3UqXjh3pL3LGTHyFTipTlBatmyJLMuYTUab56JspevNMXiqvEj/80bgHeBJrJa3Obj9MKVFOM21vQaPZP2Sb1m+dDHXbd/Oz8sW07Z1qwYRyNocuhz7Pv74Y5eAFvHx8WzYuJH/i5vKgKfexWKVq0mCcAR4nK2r38SWVc/CjYOKGfRQLj7+FpZ8+jkVpUWo1BpaRV7DiKemcUULH26/qgX+Xp71uf2CZo6vpYQfXn+M0S99gEntz5RvHyVp1gvs3/owsnwTsIoW4at49E1fQlrbEkhtP1XK8fJjRLXVcl0bf07n5ri81xfLcbFqvXqThS1HCtl9sphTx0+xZNarBLVsRebuXXz9eh5ZmbdhMUv4BZoY8VQ21/YuB+L45s0TRN3Uh0nP/ZfFP3xbrf9GTZH2LnfEyFfgQnVr7xzbtmzZzOgHHwGrFZVGg9VcTMv2HzH8vxtQafIoK7qSWU9HEtZ+IXdPeJVWkdeg0mioKC3m313pPPr8FE7qKmo8z/lQ23yxY9/HH38MVAnkYrXy26IFvHBnZyYPjbKv2WxnP3IEsBv4HAjGV7uT8e+kE/tMFm+N7UbcoC52ZxNbMJi4QV2Yek8P7u4WJhSvgPj4eP7ZvIndP3/Fta388Q8KISC4CFm+FYVyMmAgP2sgn8ddz/Z1/sj22Ct6k4VNGfl8s+kYz01+3eW9bijHxaplHPXGTZrMjbf2YfZv29h5oghdfi4zxw8mc3cqaclWYAcnDsZiMUtIigTi5hzm2t42L24JiVlzv2HlwgRu69WzRmfZqukHhUOtDUmWGyf6zqRJk+Tp06fXWiY3N9ftQu+JNrlSdb7WgdIjiBsHHWTjSkdK2EzgJSDprLIeKjXf/rkPf40HAV6etNZ6EeavdoashJpHtjWd35Fsu7r3WempotvNA5yJ0B3LIHb9vcY+6u0LxGMLmgBwGHiZXoM9uPfpV+H/27v3+Kbq+4/jryRtkkKhUGgp12K5TCaiwEEBb/CwE5XbpqDbT9Q9vFTdkLmf8EOnjjlxytAfCnMI61QcMkDnxImC4vxtqPjTowiKlx9SRO4tlNILpGnT/P7IhVJ6SUlyEvT9fDz6oGkOySfJ+eRzzvd2CIzqfuVPcxu9mHprzkgi+exmzZrFww8/bGt2owSKJJfhu5M7Te2TTpeL7xkXkN6hM8Mvn8yzDxRycN9vgMAo6D6Dqrjkp9vJG2A/oQWmKY2Ni5g2bRqFhYXcdNNNTa6rHtoGoK7uxGZum82OIyUlGMNZwKPAxcF7v6Df4Kf58fTLw2MbMts4ufj0LHIy3Cc8VjycqvtSc7msM19plYb9takuN4NHj+WOhUv50c/2cdvcr+mS6wFOA54H/g0YANjsdgaOzOfup9eyu+won++r4L3tpfxt424K39nBq5/u5/O9FRz1+po84m+qv3jy5MkA9OnTB7c7GJvTxeDRY/nVM68HrzJ0bNqQKy2dyXeYtG3/KfBPAoV3Hympv6TXgJ8wYmwqlWUHw8+bkZlNt6yO+HS1LWmgqX3y/778krUvv4j5xt+Z9/MrObjvNQKrPt0EHGDb5rY8Of0MVs7ryrX3rKJtRmZwqdLAAWrH7G6kBm+7GxkX0dKUOYCMjOO3aazwAvj9ddTW9AOWAh8RKLylwC+As+jcbQftM7Ow2WwM7dWBq43ulhXebysVX2mVhv21tTVeOrTPoF3HzgDkDTzCL58o4oqpe0lxlgMXELhO+z/w153Lji82Nfq41bU+tpVUMqx/dzq2b9vkF0pWdhdS3G3weKpxOl0cPXqU5cuXs3TpUvx+P9u2bcPjCTRt19YcG4kc6of+2dy/0vfsx/j4Xw/w17lnUlV+BlCK3fEboB/DLtnL7f/9DPn/cStV5WWUl5aQ3c7F5KHdSfFWRNRUnqj5y5IYzY1hyGzr5NMtWzhvzITgNCM/qa7nOPP8qQzN3w3AB6935Jn786k6PIdabzdSnC58NV6c7jRqa7ykOF1UezzsqoL1u32s2bKfd7Yd5MW33mfcDycdV/QvnziJp159l+Uf7GLRv7fzX0+tYfCosS28gmHYbKsIdLlcQ2CVqnlAXwacs4kRY6+g4tBBHHYb48/MYWSfTqQ4VDqipXdQWq1+f23BzTfTzl/F2d3SsdsCrSsOB4wYe4h+Z1+Lq81CbHYPMA54l4rSlayc97/4fI0/duhiBaH5kKEz65lPreWZd3fw5Pqv2bR1J8Mvv4qfz1tG5+69AbDbA+vF2ux2OnXL5ebfFTL0ksCXBsD4gj+Slr6Ap2aNZ+vGaXiq+tGuYw05uYs4Z8yd/GL+IEaMHRfeft2yJ/l6y4dsXbOEq4Z2p0t7V8R9V1pQ4LunuTEMuT17MKBn9nEtL+kZ8OM7D2O3DwSW4PfbCJwRf0Wt90X8/h9SVXGE4ZdfHR5VX3bwAPvKPWwtruSjb8rYUubggDeweEdoZsGhWge7q12UVFbj9dWFr48OgablMFsm7vSZBA6M38fvnwBUA38A+gP/SZfcTBwpqVwx9T5unLWA8YO6kqsR+zGj0c7Sao2tyVpcXIzRJoM3vyihuCKwZu0N9z8EwF3j++Cruw2YCoziS3MUd43bi82+jDsXTqRLr0BfV3lpCUsfmk5ml27HfVG526ST3rEzFdW1AFx/3+PcPWEwG1YvD8dRVxeo5v66OvoPHk7/wSNonzWEoo9zeHJme4o+aRP8goOc3h5GjjuE8YMyUp0XEDg7hyum3sfdE4cw47Izwo+7/C9Ps/wvT0c0BznR85clcVpap7ikpISCggIuHH81Cxf9ibKDgdXj7lj4CP987gk+eWcutTXTCawRfTlwOTZqcKSUU3GokvEFA064OlZ5aQmb17/OkIsncOGPruO9V58/YU3x0DZD8ydSVe7mi/ddwDjwj8NTGThjdqV56ZD9CsW7pjFktMGwS2ez6a1XqDh0gOvve5xUh53xg3Lo3iEtlm/Zd56Kr8RM53QXk4d0Z+Ouw7xXVEpdcPDTr5Ys45U/zeWTdxdQ670OuBnoj7/uTh65BbJ7VtN/SCUH9rzM9k8/p6xkb3iqUmNfKBA4Q64/AMpmd9D37NHACL7a1I8/zshl+5Y2ECy4jpQ6Bp1fzoixpfQ+4yi2RoZA9Mxsw/sffcKcB+47qUUBmlqIXgsKSP3iPO6ic9m46zCbd5XTrmNnXG3a4qv9jBTnrdR6Z5A74FGOVl5N8U4Xb6/qxNurOpHqqiN3wCEO7H6Ri38yktO+72L9qsUnTHsD8PvhaKWd/TtcvLZkK0cqHuTz9y/jSPlpx8XUNuNDJt7ahYEjK/jH4pUUf7MHp8tNTu/+9Ak+Vla6i4v6d6ar+ndjTsVXYmr//n3MvOFaHn/yKT4+aKP0iDfc9OWrKSHFuYBa7yPACLJ7PkD5wVEU73RRvNMF3AVM59D+r9iwegvvvfoZk+94BHfbOr7a5MPdpg5frQ2vx061J53DB8dSUz0Sm70//ro+fLXxTPz+wC5dsitQcL83tJIzzy/njOEVpKU3PtgkIy2V8/t2Iq9zoHnuZOcga0EBiYQr1cHw0zI5u0cG//PJDl48XNrgYPMvXHtPLn/+9bPkDZzF9k+z2FOUxlcfdwJu5m/hAc3PAYvZsLqcDavLASdt2vXkaJUDf13o6PImAI6UAxzBZvsXE245m+8PrySzSxp3T+jDsjnHRllvWL2CDatX4HS52LP/IOlulYh40TsrMRXq71w0/xEenfcY72wr5ZPdh6ksKwWbjVpv6DJqGyjemQ+kYHecR1aPWyn+pi9+/2DgdOB0/H5YOa+5Z7sNAH+wpvr9Ppzuzxg8uj3fG+qnW78SOmW7KS8t4en7pzPl7kfCUyUAUh12jNyODO6ZcdyKVNEsCqAFBSRS7lQHQ3qks37NKjbvLmfjzsN0yzsdgBf/8Fv2bV9J7wFVFO/8O5BJYErcwHo/vYG2wZ+uABwJTMHF6a7F5/scX82HwCc4Uj9n4Ag3E275Be0zD4VjCLUgbdnwJt7gFLoxY8bw2GOPqfDGWdTvrmEY+cBM0zR1LbXvsH79+jV5fdKi3cVc9OxzbN66g4Vzf8vmd9YdN9/2oitv4M/33YLffwBwEViyMvQF0xWbLZPTBl6Mp8qOIxWc7jqcrjrcbX106lpD5+5esrp5ee+132K+sQS7/SrOPO/XVFUF+oHXLVvI11sCl0C7YuqvsWGjf5d0RvbJJN11YgpEc53RU/0apcpn67lSHQzr3ZGzemTQuVNHvPXyaMPq0P60jxTni9R6/4rN7sAfHOPgSM3EV+PmrAuv5tLrbyOtbR0PTBmC13PkuOfw1cDm9Xam3D07/DcbNs7q35ttvbLZ9O9jU+jS09PVWmOBqIuvaZrrDMOYGcm2xcXFzd6fjINSFFNkXn31VebPn8/atWvDl0m79NJLuffee6mtKiPHCTln5PBW1ww+8laTkhq4XJ8j1c3bLy+j4tABMrv2IjOnO9s2vY+/bjMpThcDzh3NJddNo13HzwCoOHSAF+bdww+nPhie3jT7J+cft0BBqOmsodDfXS4XW7du5cjhUo6csFV8JeNnV1+k+dxSLkNyvtZki6lhPO+8/TazZ89mzdq1VHs8pDpdtOuUTeneneFWo1DhBfDVlGKz2an2bCQtI/BYdyx8gbVLHmfLO2/g99fhSEklIyuHzJyeVFVVkWK306eTm4Fd08hw2yks3ss111zDNddcw3PPPceePXsa/Xz379/P1KlTeeKJJyxf8CLZPjeIPiZL2xUi+cCSbRUTUEyRys7OxusNHEF7vV6ys7MZOPD469xWVlZSUFDAj6f8lIsvOh/z9b+F7yvd+w2le78BCM91TG+fQU6P3PA2a59+lG8+/5h3X1rCFVMDq081HHwVOqOu9nj4bMObOBwOfD7fcQOgEvn+JeNn11qRvoZkfK3JFlP9eLKzs8nOzqbGGziYrPFWU7p3Z6P/L7Sfj7tpxnHdKW3btiW9fQbgD+fRoHMu4J4H59I7M41uHdKO62Z56aWXwr+PGjWqyZWbZs+ezQcffMDixYubXEkrnpLtc4PoYlKjvsRMJP2d9Ztlt237itt/OYM31rwSOMp3uWnTLoO+Zw8PT50o3b+bP864np1fbm707DYl1clDL28kvX17ar3VOF0uvNUeNr61OrytLzip+OjRoxoAJUkvlEcTJ05k2rRpbN++PXzw2K1bN4qKinA6nXi9Xvp2z+LK887AbrOxe89eZt5+MzMefgL/kcNMuPo6fnrDDbyycikHSoq5oG+nk4pHU+jio8XiaxhGQSN/LjJNc10c4pFTWGv7O7t27UrXrI7UBM+Wq71eLrn0MmbcPwdfnZ8r80cyd9ZdvLziWS6ZMAlfbS3r33yNao8Hp8vNkAsuYfq99zPs+z1423aEgoICbrzxRubPn89bb73FwYMH8Xg8OBwO8vPzycrKYv/+/XG7UsypQPmc/Orn0ejRoykqKgr3x/p8vvB+HjrAPT2nHQB/eHABn370PutXLuJfa1aFH+OyC8894TlaQ1Po4qPF4mua5uKWthE5WY2dLQ/qnnHC0fbaVc+Hfw81aw/Ky2Hc8AHA8V9YhYWF3H777RQWFoa3zc3NDTeVTZs2LbwCVSKazxJJ+XxqaSw/Qge2oX+bOjN1uVwcPnw46hg0hS4+ol5e0jCMSYF/jEkxiEe+Y5pasrGxxeq7d+/OlClTIroMYXFx8QnbRrIQ/Xed8tkaofW/Wxq4FsmSpg1zxeEILLUauthILER7CVA5USxGO79AY9eNE4lCw6Ntj8eDz+dj9uzZ5OTktNisvWLFivDAkdC2e/fuVfNZC5TP1qh/zenQpf5OVihXjh4NXFAkNMZh6dKlLF26NCZ9s6f6FLpkpAsrSNKqf7Q9YMAA9u3bF9XFCtR8JonWsPUlVByjbX0JtfSMGTMmfOab1shlCCV5qPhK0lqxYgVLlizhnHPO4bPPAvN8o20qVvOZJFLDJmK32x2TArlixQoKCwvp1asXfr9fB5enAE01kqQW65GWaj6TRIp364uWNz11qPhKUlNTsXzb1C+QCxYsiGnriw4uTx0qvpL0dDQv3yb1C+SDDz6YlCs3Sfyp+ErS09G8iHzbaMCViIiIxVR8RURELKbiKyIiYjEVXxEREYup+IqIiFhMxVdERMRiKr4iIiIWU/EVERGxmIqviIiIxVR8RURELKbiKyIiYjEVXxEREYup+IqIiFgsqqsaGYbRAcgP3hxmmubM6EMSkURQPotYJ9oz36uATNM0XwAwDKMg+pBEJEGUzyIWierM1zTNxfVu5gGLmtu+uLi42ccrKyuLJpy4UEyRUUyRScaYQlqTzy3lMiTna022mJItHlBMkYo2pqiKb4hhGHlAqWmaRc1tl52d3eJjRbKN1RRTZBRTZJIxpvoiyedIX0MyvtZkiynZ4gHFFKloYmqx+DbR9FRkmua6ercnmaZ5y0lHISKWUD6LJIcWi2+DpqgTGIYxyTTN3wd/z2+QxCKSRJTPIskh2tHO+cAcwzDuDv5JoyNFTlHKZxHrRDvgah3QJ0axiEgCKZ9FrKNFNkRERCym4isiImIxFV8RERGLqfiKiIhYTMVXRETEYiq+IiIiFlPxFRERsZiKr4iIiMVUfEVERCym4isiImIxFV8RERGLqfiKiIhYTMVXRETEYiq+IiIiFlPxFRERsZiKr4iIiMVUfEVERCym4isiImIxFV8RERGLqfiKiIhYLCXaBzAMIz/46w9M05wZ7eOJSOIon0WsEVXxNQxjCMEkNQxjpmEYeaZpFjW1/axZs6J5OhGJo9bks3JZJDo2v98f9YMYhtEBmGOa5i3RhyQiiaR8Fom/WPX5GkBZjB5LRBJL+SwSZy2e+RqGUdDIn4tM01zXYLtFwBumab4Qw/hEJIaUzyLJocU+X9M0Fzd1n2EYc4BtwW3KgMwYxiYiMaZ8FkkOUfX5GoaRB+QFb05WH5HIqUv5LGKdmAy4EhERkchpkQ0RERGLRb3IxskyDGMSgX6lIaZp/r6191sdU3D6RWgBgmFWLUAQ6ftgGMYcK2KK4HMbQrDp0qrBOq3Yl/Ka6/OMcUz5wEzTNH/QxP2W79/xolyOPqYG21mSy5HEpHwOP2fM8zkhZ77BQAmOsCyrt6pORPcnIibgKiAztAM2MWrU6phC2+VzrK8u0fHcEnyP8oJ9iAmNKXg7NJq3KPhlEncNRw83iMny/TtelMsxiym0nSW53IqYlM/EJ58T1ew8DAitnFMENHwDW7rf8phM01xc7ygrD2jyw7AqJggPkmlyVTEr4wl+iX0YXBnp982tdmZVTIAJPB86gjdN8yMLYmpJIvbveFEuxyAmsDyXW4xJ+Ryxk9rHE1V8OzS43amV98dDRM8ZTJBSi3bESGJqdknPGGspnj7Bn1LDMBYFm/cSGpNpmmXAIuD5YGzJIBH7d7wolyOTbLkMyudYOal9PFHFt6U5hImYYxjpc06ycApGszEZhpHfXHOI1fEEbQsmyIdA3JvzaPk9mgSsM02zT73bifZtmkOrXI5MsuUyKJ9j5aT28UQV3w84drSQB7zRyvsTEROGYUwKdaZb1E/XUkylhmHkB3fAPAv6PyL53EI6YM0ShS3FVL9p6iGSo+glYv+OF+VybGKyOpcjiUn5HJmT2scTUnzrdeDnB2+vAzAM443m7k9kTMG/zzEM40PDMD6MdzyRxGSa5kfBv2VyYtNHIuJ5AehQ7/64j0RsKSZgsWEYBcH7r7JwdOSkwD/HjswTuX/Hi3I5NjFZncsRxqR8DopHPmuRDREREYtpkQ0RERGLqfiKiIhYTMVXRETEYiq+IiIiFlPxFRERsZiKr4iIiMVUfEVERCz2/xuTcvfPDaZtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set into eval mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Initialize plots\n",
    "f, (y1_ax, y2_ax) = plt.subplots(1, 2, figsize=(8, 3))\n",
    "\n",
    "# Make predictions\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    test_x = torch.linspace(0, 1, 51)\n",
    "    predictions = likelihood(model(test_x))\n",
    "    mean = predictions.mean\n",
    "    lower, upper = predictions.confidence_region()\n",
    "    \n",
    "# This contains predictions for both tasks, flattened out\n",
    "# The first half of the predictions is for the first task\n",
    "# The second half is for the second task\n",
    "\n",
    "# Plot training data as black stars\n",
    "y1_ax.plot(train_x.detach().numpy(), train_y[:, 0].detach().numpy(), 'k*')\n",
    "# Predictive mean as blue line\n",
    "y1_ax.plot(test_x.numpy(), mean[:, 0].numpy(), 'b')\n",
    "# Shade in confidence \n",
    "y1_ax.fill_between(test_x.numpy(), lower[:, 0].numpy(), upper[:, 0].numpy(), alpha=0.5)\n",
    "y1_ax.set_ylim([-3, 3])\n",
    "y1_ax.legend(['Observed Data', 'Mean', 'Confidence'])\n",
    "y1_ax.set_title('Observed Values (Likelihood)')\n",
    "\n",
    "# Plot training data as black stars\n",
    "y2_ax.plot(train_x.detach().numpy(), train_y[:, 1].detach().numpy(), 'k*')\n",
    "# Predictive mean as blue line\n",
    "y2_ax.plot(test_x.numpy(), mean[:, 1].numpy(), 'b')\n",
    "# Shade in confidence \n",
    "y2_ax.fill_between(test_x.numpy(), lower[:, 1].numpy(), upper[:, 1].numpy(), alpha=0.5)\n",
    "y2_ax.set_ylim([-3, 3])\n",
    "y2_ax.legend(['Observed Data', 'Mean', 'Confidence'])\n",
    "y2_ax.set_title('Observed Values (Likelihood)')\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
